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Abstract 

We present lattice-gas modeling of the steady-state behavior in CO oxidation 
on the facets of nanoscale metal clusters, with coupling via inter-facet CO dif- 
fusion. The model incorporates the key aspects of reaction process, such as 
rapid CO mobility within each facet, and strong nearest-neighbor repulsion 
between adsorbed O. The former justifies our use a "hybrid" simulation ap- 
proach treating the CO coverage as a mean-field parameter. For an isolated 
facet, there is one bistable region where the system can exist in either a re- 
active state (with high oxygen coverage) or a (nearly CO-poisoned) inactive 
state. Diffusion between two facets is shown to induce complex multistabil- 
ity in the steady states of the system. The bifurcation diagram exhibits two 
regions with bistabilities due to the difference between adsorption properties 
of the facets. We explore the role of enhanced fluctuations in the proximity 
of a cusp bifurcation point associated with one facet in producing transitions 
between stable states on that facet, as well as their influence on fluctuations 
on the other facet. The results are expected to shed more light on the reaction 



1 



kinetics for supported catalysts. 
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I. INTRODUCTION 



A central problem in the modeling of catalytic surface reactions is to bridge the gap 
between the studies of reactions on ideal homogeneous surfaces of macroscopic size and 
real catalytic processes on supported catalysts. The latter are composed of small nm-sized 
metallic clusters on an inert support material. The surface of these clusters contains facets 
of different orientations, and hence differing reactivities, which are coupled, e.g., through 
diffusion of mobile adsorbates. Due to this coupling effect, and since fluctuations in such 
small size facets can become dominant, the behavior of these systems can differ strongly 
from that predicted by mean-field rate laws on ideal single- crystal surface planes. An un- 
derstanding of the effects that can arise is crucial not only for the modeling of real catalysis, 
but also for the design of nanostructured catalysts. 

For catalytic CO oxidation, by far the majority of experimental studies have been per- 
formed on extended single-crystal substrates. Mean-field equations, ignoring spatial corre- 
lations and ordering in surface adlayers, have been successfully applied to predict a wide 
variety of surface phenomena, including nonlinear dynamics and spatial pattern formation 
imi. More sophisticated lattice-gas models can be utilized to describe the effect of adspecies 
interactions, and "hybrid" formulations can treat directly and efficiently rapid surface diffu- 
sion of some adsorbates such as CO (the rate of CO diffusion is of many orders of magnitude 
higher than the rates of adsorption, desorption and O diffusion). A recently developed 
canonical hybrid lattice-gas model of this type for CO-oxidation incorporates a random 
distribution of CO, and strong nearest-neighbor O repulsion 

In-situ experimental studies of reactions on supported catalysts are still beyond the 
imaging capabilities of current microscopes. However, for a reasonable model system as far 
as the dimension of the facets and the coupling of different orientations are concerned, one 
can consider reactions on a field emitter tip. The surface of such a tip can be imaged with 
field electron/field ion microscopy with nearly nanometer/atomic resolution. In a series of 
recent experiments, the spatio-temporal dynamics of the fluctuations which arise in catalytic 
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CO oxidation on a Pt tip have been studied in some detail 0. 

Conventional or hybrid lattice-gas modeling is ideally suited to analyze both kinetics and 
fluctuations in nanoscale systems 0,0. For example, recent modeling of reaction kinetics 
on the facets coupled by reactant diffusion reveals shifts in the reaction window for bistable 
systems, and shifts in the oscillation window for systems which incorporate additional feed- 
back. The latter also indicates a possible transition from regular to chaotic oscillations due 
to additional supply of reactants from adjacent facet Lattice-gas modeling has also 
been applied to analyze fluctuations in catalytic CO oxidation within a small isolated facet 
containing hundreds to thousands of adsorption sites What has been neglected in these 
studies is the coupling between adjacent facets on the Pt field-emitter tip. Experimen- 
tally, there is some justification for this, since one observes that fluctuations on different 
facets were usually largely uncorrelated (while a high degree of spatial correlations exists 
only for fluctuations on a single facet). One can in fact explain this behavior noting that 
CO diffusion is very fast within a facet but inhibited by structural heterogeneities (such 
as atomic steps) at the periphery of the facet. However, there is some coupling between 
facets, as the kinetics is described by a global bifurcation diagram, and there must be some 
correlation between the fluctuations. 

Thus, to shed more light on reaction behavior for coupled nanoscale-facets, we focus in 
this work on incorporation of inter-facet diffusion into the above-mentioned canonical model 
describing key aspects of CO oxidation. We find that the diffusive coupling between two 
adjacent facets with different adsorption properties results in the appearance of complex 
multistability of the steady states (which is fully characterized exploiting a novel simulation 
approach). Consequently, the bifurcation diagram exhibits two regions with bistabilities 
instead of single bistable region observed in the case of isolated facet. We show that typically 
the different facets have different fluctuation characteristics. However, in some cases, the 
diffusion-induced coupling can change radically the fluctuation behavior of whole system, 
producing a transition between stable states on one facet induced by enhanced fluctuations 
on the other facet. 
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II. THE TWO-FACET REACTION MODEL AND ITS ANALYSIS 



A. Model Specification 

The lattice-gas model employed for our studies is based on Langmuir-Hinshelwood mech- 
anism for catalytic CO oxidation. The two "adjacent" small facets labeled i = 1,2 repre- 
sented by two rectangular grids of Ni = Lq x Lj adsorption sites. One could physically 
connect the facets, e.g., at a common edge of length Lq, which would introduce spatial 
inhomogeneity within each facet. However, we adopt a simplified treatment invoking peri- 
odic boundary conditions. We incorporate the following steps into a "hybrid" model [0,0, 
adapted here for the two-facet case: 

(i) CO (gas) adsorbs onto single empty sites at rate pco per site on both facets, and desorbs 
from at rate d. CO (ads) can hop very rapidly to other empty sites within each facet. Since 
we consider below the regime of infinitely mobile CO (ads) within the facets, and neglect 
interactions between CO (ads) and other CO (ads) and O(ads), the distribution of CO (ads) 
on sites not occupied by O(ads) is random. Thus, in simulation procedure, we track only 
the total numbers N^q and N^q of CO (ads) on both facets, whereas distribution of O(ads) 
for each facet is described by a full lattice-gas approach. 

(ii) CO (ads) can diffuse at a finite rate from each facet to empty sites on the adjacent facet. 
Physically, CO (ads) would hop across a common edge (of length Lq) to adjacent empty sites 
at rate, say, h'. In our model with periodic boundary conditions, this mass transport is 
mimicked by transport of CO (ads) from facet i to j at rate (in molecules per unit time) 
hiN^Q times the probability that a site is empty on facet j, where hi = h' /Li. 

(iii) 02(gas) adsorbs dissociatively at the impingement rate po^ per site (as described below), 
but the sticking coefficients for oxygen differ for each facet: Sq 7^ Sq so that the two facets ex- 
hibit a different catalytic activity. To account for very strong NN 0(ads)-0(ads) repulsions, 
we invoke an "eight-site rule" wherein 02(gas) adsorbs only at diagonal nearest-neighbor 
(NN) empty sites, provided that the additional six sites adjacent to these are not occupied 
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by O(ads). Also, O(ads) is immobile, and can not desorb, and thus never occupies adjacent 
sites within each facet. 

(iv) Adjacent pairs of CO (ads) and O(ads) within a facet reacts at rate k to form C02(gas). 

As is the single-facet case, we set pco + P02 = k = 1 and consider system behavior as a 
function of pco with = 1 — Pco- In contrast to the simplified ZGB (Ziff-Gulari-Barshad) 
model PI, where CO(ads) are immobile, the main feature of the above model when applied 
to single isolated facet is the appearance of strong bistability: over a significant range of 
Pco and (smaller) d values, there are both stable reactive state (with high O coverage 0o) 
and stable inactive state (with the surface nearly poisoned by CO). 

B. (Standard) Constant-Partial-Pressure and (Refined) Constant- Coverage 

Simulations 

Standard "constant (partial) pressure" kinetic Monte Carlo (KMC) simulations simply 
implement all processes (adsorption, desorption, reaction, diffusion) with the appropriate 
relative rates. They can be readily implemented for single and coupled facet models, and 
enable analysis of stable steady states as well as kinetics. To also analyze the non-trivial 
unstable steady states existing in bistable (or multistable) regions, a constant-coverage algo- 
rithm has been applied previously for large single facet systems. In the adsorption part 
of the algorithm, instead of attempting to deposit with the appropriate relative rates, these 
simulations specify a target CO coverage Oco, and deposit CO (O) for actual CO coverages 
below (above) this target. Then the asymptotic pco value (corresponding to given 9co) is 
found as a fraction of attempts to deposit CO. In our case of two communicating facets with 
different O sticking properties, the above const ant- coverage approach is unsuitable: fixing 
the CO coverage on both facets corresponds to different values of pco (which is unphys- 
ical). To resolve this problem, in this work, we develop a combined constant-pressure & 
constant-coverage algorithm. Analysis using this novel simulation procedure is performed 
in two steps. 
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First, we fix a target CO coverage for tlie first facet Qqo, but use tlie constant-pressure 
algoritlim for the second facet. Fixing Qqo during this step allows us to probe both stable 
states together with unstable state in the bistability region on the first facet, but the CO 
coverage on the second facet Qqq relaxes with time to a stable steady-state value corre- 
sponding to the rate pco (which is equal for both facets). Because of this feature, we can 
not probe the unstable state on the second facet in the bistable region (the same limitation 
as in constant-pressure simulations). 

Second, to probe the full set of states on the second facet (including the unstable state), 
we invert the above procedure, prescribing the coverage Qqq and performing standard 
constant-pressure simulations on the first facet. This gives us all states for the second facet, 
but just the stable states for the first one (matching results from above). Thus, whereas the 
first step gives us the full picture for the first facet and incomplete data for the second one, 
after performing both steps, we can obtain the full picture (i.e., the complete bifurcation 
diagram) for the whole system. 



C. Analytic Pair- Approximation Treatment 

Our novel simulation procedure extends conventional KMC approaches in order to obtain 
complete steady state behavior. It is thus natural to compare results with those obtained 
from analysis of (approximate) rate equations for our two-facet CO oxidation (which auto- 
matically provide complete steady state information). For convenience, we denote adsorbed 
CO by A and adsorbed O by B, using standard notation for monomer-dimer A + B2 reaction 
model. Following 0, we let E denote empty sites, Z = 1 — B denote sites not occupied by 
B, and indicate probabilities of various configurations by square parentheses. Then, for CO 
and O coverages on the ith facet (z = 1,2) of Qqq = [Ai] and Q^o ~ [-^«]' respectively, the 
rate equations have the following form: 

^\Ai] = pco[E^] - d[Ai] - Ak[AiB,] - h,Ji, (1) 
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The net diffusive CO flux Ji = -J2 = [A1E2] - [A2E1] = - ^2)] - [^2(1 - Si)] 

introduced in (|ip describes the CO transport between facets with the microscopic hop rate 
h' related to the hi as described above. (In our analysis, we set A^^i = N2 corresponding to 
equal size facets, and thus hi = h2 = h.) The lack of spatial correlations between facets 
further reduces the J, to Ji = -J2 = [Ai]{l - [B2]) - [A2](l - [5i]). The lack of NN O(ads), 
together with the random distribution of CO within each facet, allows further reduction of 
(0). For example, the reaction terms can be exactly expressed in terms of the coverages 0: 
[AB,] = [Ai][Bi]/{l-[B,]). 

These results allow derivation of an exact relationship between [A] and [B] for the steady 
states {d/dt[Ai] = 0) from the first equation in (|T]): 



[Ai] = -pco(l - [5i])(/^i(l - [^2]) + h)/c 
[A2] = -Pco(l - [B2]){h2{l - m) + b2)/c, 



(2) 



where c = hih2{l - [Bi]){l - [^2]) -6162; h = Pco + d + 4k[B2]/{l - [B2]) + /i2(l - [5i]) and 
&2 = Pco+'^+4/i;[-Bi]/(l — [i?i]) + /;.i(l — [-B2])- We next substitute the resulting expressions (0) 
for [Ai] into the second equation (|ip and applying the pair (Kirkwood-type) approximation 
to decouple the eight-site correlation functions in (||) describing the "eight-site rule". 
The analysis of the steady states for two adjacent facets (corresponding to the conditions 
d/dt[Ai] = and d/dt[Bi] = 0) is then reduced to the solving the closed system of two high- 
order polynomial equations with respect to [Bi], and subsequent evaluation of [Ai] from 



III. RESULTS FOR STEADY-STATE BIFURCATION BEHAVIOR 
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A. Case Study: sj) = 0.2 and = 1.0. 



We now present a brief overview of new results from our novel KMC simulation procedure 
for coupled Na x N2 site facets, which are compared with pair-approximation predictions 
using rate equations (|l|). To demonstrate the effect of inter-facet diffusion, for h = 0.01, we 
show in Fig. 1 the CO coverages vs. pco for = 0.2 and Sq = 1.0. A very good agreement 
is found between results from KMC simulation and the pair-approximation (|T]). The latter 
produces only slight deviation of bistable regions towards lower pco- The difference between 
[Ai] and [Bi] obtained using both approaches is the most significant in the reactive state, 
since the high O coverage makes the pair-approximation for "eight-site-rule" less accurate. 
In contrast to the case of isolated facets when h = (shown in the insets of Fig. 1), a 
qualitatively new feature appears in Qqq and Qqq behavior due to the inter-facet communi- 
cation by CO hopping. Besides the bistability region (#1) for 0.113 < pco < 0.1515 on the 
first facet with Sq = 0.2 (resembling closely that for h = 0), an additional bistable region 
(#2) at 0.17 < Pco < 0.3614 forms on the first facet covered predominantly by CO. This 
is a response to the bistable behavior in this range of pco values occurring on the second 
facet with Sq = 1.0. Analogously, an additional bistable region on the second facet for 
0.113 < Pco < 0.1515 arises in the state with low CO coverage as a response to a bistable 
behavior on the first facet. 

Fig. 2 shows the variation with h of the bifurcation diagrams in the {pco, c?)-plane. 
For low diffusion rate h, the phase diagram contains two overlapping bistable regions (see 
Fig. 2(a)): the small region ^1 at lower pco terminating at cusp point {pi ~ 0.1548, ^ 
0.025) originates from the first facet {sq = 0.2), whereas the larger region #2 terminating at 
cusp point {pI ~ 0.409, ~ 0.051) occurs due to the second facet. Note that the behavior 
of the both bistable regions differs dramatically with varying h. The larger region #2 varies 
only weakly with h, whereas the smaller region 7^1 first shrinks, and then disappears. The 
latter behavior can be understood as CO on the first facet is supplied very rapidly to the 
second facet with lower CO coverage, thus suppressing completely the occurrence of region 
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B. Dependence of Behavior on Sticking Probabilities 

In extending our analysis to general sq values, we introduce three different coupling 
regimes: {i) < h < hi, weak-coupling regime, when two overlapping bistable regions cross 
each other on the diagram {pco, d); (ii) hi < h < h2, medium-coupling regime with one 
smaller bistability locating completely inside another larger one; (iii) h > h2, strong- coupling 
regime with only one common bistable region. This classification will be valuable for the 
studies of fiuctuation behavior below. The values hi and /i2 separating these regimes are 
shown in Fig. 3(a) as the functions of Sq (at Sq = 1.0). As Sq — > 0, obviously both hi and 
/i2 — (since the first facet becomes an oxygen- free "sink" for CO, and thus the bistabihty 
region #1 disappears, see the inset of Fig. 3(a)). 

The variation of hi and /i2 with Sq shows three distinct branches arising from competi- 
tion between the effects of: (i) the inequality in oxygen sticking, and (ii) finite inter-facet 
transport. The first branch (sq < 0.65) is formed due to the increased Sq value providing 
the expansion of the first smaller bistable region #1. As the difference in oxygen sticking 
properties jsg — SqI decreases, the second branch (0.65 < Sq < 1) appears due to a dominant 
role of restricted inter- facet transport. We comment in particular on the case Sq = Sq — 1 
(identical facets) where h2 ~ 0.06 ^ 0, and where one might expect just one bistable region 
(so /i2 = 0). Indeed, the "symmetric" steady-states for a single facet must apply to this 
system, but restricted inter-facet transport produces two additional " unsymmetric" stable 
steady-states (see Fig. 3(b), inset). The region of the existence of these two states located 
inside the large bistable region #2, is indicated on the bifurcation diagram in Fig. 3(b). For 
the third branch {sq > 1), the inequality in sticking again begins to be of primary importance 
producing the increase in hi and /i2 shown in Fig. 3(a). 
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IV. RESULTS FOR FLUCTUATION BEHAVIOR 



Bearing in mind the classification of coupling regimes in Sec.IIIB, we turn now to the 
analysis of the coverage fluctuations. The fluctuation amplitudes depend inversely on system 
size, and become more significant approaching a cusp critical point (analogous to fluctuation 
behavior near a critical point in equilibrium thermodynamics). Thus, they can induce 
transitions between the stable steady states for systems which are small enough, or in a 
close proximity to a cusp point [§,|T3. Here, we consider two coupled facets, each with 
30x30 lattice sites, and analyze fluctuation effects both for different coupling regimes, and 
for different Sq values at fixed Sq = 1.0. As already noted, the size and location of the small 
bistability region #1 depends strongly on the value of |sq — Sq]. We show below that this 
factor, together with the coupling strength {h), strongly influences fluctuation behavior. We 
consider separately two regimes of coupling strength. 



A. Weak Coupling Regime 

In Fig. 4, we plot the normalized amplitudes F] = [Ni{5{Qiy) /{{&]){!- {Q])))]^^'^ 
(j = CO, O) ((e;) denote the tome averages of = N'^/Ni and 50) = B) - {0])) vs. 
d on both facets close to cusp point {pl,dl) terminating the large bistability region #2 for 
Sq = 0.7 and h = 0.005 (in this case ~ 0.4096). As d approaches d^, the amplitudes 
increase drastically, whereas just slight increase of Fj is observed with d ^ d"^, reflecting 
very weak response of the first facet on the critical behavior of the second facet. Moreover, 
we find that correlations between fluctuations on the different facets (measured by Fj-"^ = 
(6]6|) — (0])(0^) - Fig. 4, insets) are practically zero compared to correlations within each 
facet (measured by Fj). An analogous (but reversed) picture emerges close to the other cusp 
point. This behavior is consistent with the observed weak correlation between fluctuations 
on different facets of Pt field-emitter tips . 

The distinct character of fluctuations on different facets is also clear when analyzing time 
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series for ©j and or the normalized probability distributions P(©^) (with J dxP{x) — 1). 
Fig. 5 shows P(6* ) obtained near the midpoints of the first and second bistable regions close 
to cusp points. As d ^ dl (Fig. 5a), the distributions of Qqq and Qq become bimodal, 
reflecting a transition between reactive and inactive states on the first facet induced by 
increased fluctuations. However, the form of Qqq and ©q distributions remains effectively 
mono- modal. Since the difference in ©| in the two stable states of bistability region #1 on 
the second facet is very small (see coverages vs. pco in Fig- 1), the transition between these 
states is practically invisible in the ©| time series (Fig. 5(a)). Analogous (but reversed) 
behavior is observed as d ^ d^ (Fig. 5(b)). 



B. Medium Coupling Regime 

To study what happens when the cusp point {pi, d]) nears {pi, d^.) in the {pco, (i)-plane, 
we consider next the medium-coupling regime when the difference |so — SqI i^ small. We set 
Sq — 0.9, h — 0.018 , and analyze the ©* time series close to the first cusp point (pj 0.383, 
dl ~ 0.04) of the bifurcation diagram shown in Fig. 6(a). Specifically, we examine the 
midpoint of the larger bistability region #2, close to the upper boundary of the inner region 
#1 (Fig. 6(b), d — 0.035 and pco ~ 0.375). Besides the conventional fluctuation-induced 
transition between two states with low and high Qqq in bistability region #1 (for which Qqq 
is always low, and thus distinct transitions in this quantity are masked), the time series for 
©CO show the occurrence of additional transition between two states with high ©co- These 
derive from bistability #2, as is clear from the corresponding transitions in the time series for 
©^. This produces a rather complicated form of the distributions for ©co) with additional 
peak splitting. Since the second cusp point (p^, d^) is close, corresponding fluctuations of 
the coverages in the states of region #2 on the first facet are sufficiently large to induce 
the additional transition. Note that the time series shown on Fig. 6(b) are synchronized for 
both facets due to the inter-facet CO diffusion. 
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V. SUMMARY 



We have shown that the diffusion of adsorbed CO between two small adjacent facets 
in catalytic CO oxidation leads to the appearance of complex multist ability on each facet. 
This replaces the conventional simpler bistability with coexisting reactive and inactive stable 
steady states. The additional stable steady states arise due to the difference in the sticking 
properties of the facets, and the finite rate of inter-facet CO diffusion in our model (whereas 
the CO mobility within each facet is assumed to be infinite). The topology of bifurcation 
diagrams (pcO) d), ^ well as the effect of coverage fluctuations, are shown to be qualitatively 
different in regimes with different strengths of inter-facet coupling. In the weak-coupling 
regime the correlations between the fluctuations near the cusp points terminating bistable 
regions are low. However, in the medium-coupling regime, the transitions between the 
additional bistable states on each facet can occur due to the communication between facets 
via CO diffusion. 



ACKNOWLEDGEMENTS 



N.P. gratefully acknowledges the support of this research by Alexander von Humboldt 
Foundation and thanks the Ames Laboratory for warm hospitality and providing with com- 
putational facilities. D.-J.L. and J.W.E. were supported for this work by the Division of 
Chemical Sciences, USDOE, through Ames Laboratory (operated by for the USDOE by 
Iowa State University under Contract No. W-7405-Eng-82). 



13 



REFERENCES 

* Present address: Institut fiir Theoretische Physik T34, Physik-Department der TU 
Miinchen, James-Pranck-Strasse, D-85747 Garching b. Miinchen 

[1] R. Kapral and K. Showalter, Chemical Waves and Patterns (Kluwer, Dordrecht, 1995); 
R. Imbihl and G. Ertl, Chem. Rev. 95, 697 (1995). 

[2] M. Tamarro, M. Sabella, and J.W. Evans, J. Chem. Phys. 103, 10277 (1995). 

[3] E.W. James, C. Song, and J.W. Evans, J. Chem. Phys. Ill, 6579 (1999). 

[4] V.P. Zhdanov and B. Kasemo, Surf. Sci. 412/413, 527 (1998). 

[5] Yu. Suchorski, J. Beben, E.W. James, J.W. Evans, and R. Imbihl, Phys. Rev. Lett. 
82, 1907 (1999); Yu. Suchorski, J. Beben, R. Imbihl, E.W. James, Da- Jiang Liu, and 
J.W. Evans, Phys. Rev. B 63, 165417 (2001). 

[6] V.P. Zhdanov and B. Kasemo, Surf. Sci. Rep. 39, 25 (2000). 

[7] C.R. Brundle, R.J. Behm, and J. A. Barker, J. Vac. Sci. Technol. A 2, 1038 (1984); 
S.-L. Chang and P.A. Thiel, Phys. Rev. Lett. 59, 296 (1987). 

[8] R.M. Ziff, E. Gulari, and Y. Barshad, Phys. Rev. Lett. 56, 2553 (1986). 

[9] R.M. Ziff and B.J. Brosilow, Phys. Rev. A 46, 4630 (1992). 

[10] K. Fichthorn, E. Gulari, and R. Ziff, Phys. Rev. Lett. 63, 1527 (1989). 



14 



FIGURE CAPTIONS 



Fig. 1. pco-dependences of CO coverages for the facets with (a) Sq — 0.2 and (b) Sq — 1.0 
at d = 0.005 and h = 0.01. The insets show the corresponding dependences for the isolated 
facets obtained using pair- approximation. 

Fig. 2. Bifurcation diagrams (pco, d) showing the bistable regions #1 and #2 calculated 
for different coupling regimes at Sq = 0.2 and Sq = 1.0. 

Fig. 3. Dependence of the critical parameters hi and /i2 on the Sq at fixed Sq — 1.0. The 

insets show pco-dependences of the steady-state CO coverages on the two facets at Sq = 0, 
Sq — 1.0 and h — 0.02. (b) Bifurcation diagram calculated using pair approximation at 
Sq — Sq — 1.0 and h — 0.01 (medium-coupling regime). #1 and #2 indicate bistable region 
^1 (where both " unsymmetric'' steady states are stable) and bistable region 7^2, respec- 
tively. The inset shows the corresponding pco-dependence of O^q — ©co — ®co obtained 
a,t d — 0.005 with the stable and unstable states indicated by black and gray colors respec- 
tively. 

Fig. 4. Normalized amplitudes F] = [Ar.((J(e})2)/ ((ep(l - (Op))] and correlations 
= (0]0j) — (0])(0j) of coverages fluctuations on weakly coupled facets for species 
j = CO and O vs. d for pco ~ 0.4096 close to the cusp point of bistability ^1; L = 60 and 
700000 MC steps are taken for combined MC procedure. 

Fig. 5. Time series and distributions of 0* close to midpoint of (a) bistable region 7^1, (b) 
bistable region #2 {h = 0.005, = 0.7, = 1.0). 

Fig. 6. (a) Bifurcation diagram calculated using combined MC scheme at Sq = 0.9, Sq = 1.0 
and h = 0.018 for the facets with medium coupling strength. ^1 and ^2 indicate bistable 
regions #1 and #2, respectively. The inset shows the stationary solutions ©co vs. pco 
on each facet at d = 0.035. (b) 0* distributions and corresponding time series obtained 
at d = 0.035 and pco = 0.375 close to the cusp point of bistability ^1 in the bifurcation 
diagram (a). 
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